knitr::opts_chunk$set(fig.width = 6, fig.height = 4, fig.path = 'figs/',
echo = TRUE, message = FALSE, warning = FALSE)
library(oharac) ### remotes::install_github('oharac/oharac')
oharac::setup()
source(here('common_fxns.R'))
options(dplyr.summarise.inform = FALSE)
clean_titles <- function(df) {
x <- tibble::tribble(
~raw, ~clean,
"cephalopods", "Cephalopods",
"corals", "Corals",
"crustacea_arthropods", "Crust/Arthro",
"echinoderms", "Echinoderms",
"elasmobranchs", "Elasmobranchs",
"fish", "Fish",
"marine_mammals", "Mammals",
"molluscs", "Molluscs",
"polychaetes", "Polychaetes",
"reptiles", "Reptiles",
"seabirds", "Seabirds",
"sponges", "Sponges",
"air_temp", "Air Temp",
"biomass_removal", "Biomass Removal",
"bycatch", "Bycatch",
"entanglement_macroplastic", "Macroplastic",
"eutrophication_nutrient_pollution", "Nutr Poll",
"habitat_loss_degradation", "Hab Loss Degr",
"inorganic_pollution", "Inorg Poll",
"invasive_species", "Inv Species",
"light_pollution", "Light Poll",
"noise_pollution", "Noise Poll",
"oa", "OA",
"oceanographic", "Oceanographic",
"organic_pollution", "Org Poll",
"plastic_pollution_microplastic", "Microplastic",
"poisons_toxins", "Poisons Toxins",
"salinity", "Salinity",
"sedimentation", "Sedimentation",
"slr", "Slr",
"storm_disturbance", "Storm Dist",
"uv", "Uv",
"water_temp", "Water Temp",
"wildlife_strike", "Wildlife Strike"
)
df <- df %>%
left_join(x %>% rename(tx_clean = clean), by = c('taxon' = 'raw')) %>%
left_join(x %>% rename(str_clean = clean), by = c('stressor' = 'raw'))
return(df)
}
vuln_tx_df <- read_csv(here('_output', 'vuln_gapfilled_tx.csv'))
vuln_score_df <- read_csv(here('_output', 'vuln_gapfilled_score.csv')) %>%
gather(stressor, vuln, -vuln_gf_id)
vuln_sd_df <- read_csv(here('_output', 'vuln_gapfilled_sd.csv')) %>%
gather(stressor, vuln_sd, -vuln_gf_id)
vuln_df <- vuln_tx_df %>%
full_join(vuln_score_df, by = 'vuln_gf_id') %>%
full_join(vuln_sd_df, by = c('vuln_gf_id', 'stressor')) %>%
clean_titles()
taxa <- vuln_df$taxon %>% unique() %>% sort()
strs <- vuln_df$stressor %>% unique() %>% sort()
test_df <- vuln_df %>%
filter(taxon %in% taxa[1:6]) %>%
filter(stressor %in% strs[1:6]) %>%
rowwise() %>%
mutate(vuln_d = rnorm(n = 1, mean = vuln, sd = vuln_sd)) %>%
group_by(stressor, taxon) %>%
mutate(n_inv = 1/n()) %>%
ungroup() %>%
mutate(n_inv = n_inv / max(n_inv))
mean_df <- test_df %>%
group_by(tx_clean, str_clean) %>%
summarize(mean = mean(vuln, na.rm = TRUE))
linecol <- 'grey20'; linesize = .2
x <- ggplot(test_df) +
theme_ohara() +
geom_vline(aes(xintercept = vuln_d, alpha = n_inv),
color = linecol, size = linesize) +
geom_hline(aes(yintercept = vuln_d, alpha = n_inv),
color = linecol, size = linesize) +
geom_vline(data = mean_df, aes(xintercept = mean), color = 'yellow', size = .3) +
geom_hline(data = mean_df, aes(yintercept = mean), color = 'yellow', size = .3) +
geom_vline(data = mean_df, aes(xintercept = mean), color = 'red', size = .2) +
geom_hline(data = mean_df, aes(yintercept = mean), color = 'red', size = .2) +
geom_point(data = mean_df, aes(x = mean, y = mean), color = 'red') +
coord_fixed(xlim = c(0, 1), ylim = c(0, 1), expand = TRUE) +
scale_x_continuous(breaks = c(0.1, 0.9), labels = c(0, 1)) +
scale_y_continuous(breaks = c(0.1, 0.9), labels = c(0, 1)) +
scale_alpha(range = c(0.02, 0.06)) +
facet_grid(tx_clean ~ str_clean) +
theme(axis.title = element_blank(),
legend.position = 'none',
panel.grid.major = element_blank(),
strip.background = element_blank(),
strip.text.x = element_text(angle = 90, hjust = 0),
strip.text.y = element_text(angle = 0, hjust = 0))
ggsave('vuln_summary_plot.png', width = 6, height = 6, dpi = 300)
knitr::include_graphics('vuln_summary_plot.png')

taxa <- vuln_df$taxon %>% unique() %>% sort()
strs <- vuln_df$stressor %>% unique() %>% sort()
test_df <- vuln_df %>%
filter(taxon %in% taxa[1:6]) %>%
filter(stressor %in% strs[1:6]) %>%
rowwise() %>%
mutate(vuln_d = rnorm(n = 1, mean = vuln, sd = vuln_sd)) %>%
ungroup()
zeros <- test_df %>%
group_by(tx_clean, str_clean) %>%
summarize(all_zero = all(vuln == 0), .groups = 'drop')
mean_df <- test_df %>%
group_by(tx_clean, str_clean) %>%
summarize(mean = mean(vuln, na.rm = TRUE), .groups = 'drop')
dens_df <- test_df %>%
group_by(tx_clean, str_clean) %>%
summarize(x = list(KernSmooth::bkde(x = vuln_d,
kernel = 'normal',
bandwidth = .05) %>%
as.data.frame())) %>%
unnest(x) %>%
### constrain x from 0 to 1
filter(between(x, -.02, 1.02)) %>%
### rescale y from 0 to 1 for each plot; change grouping to allow for relative scales
group_by(tx_clean, str_clean) %>%
complete(x = c(0, 1), fill = list(y = 0)) %>%
mutate(y = (y - min(y))/ (max(y) - min(y))) %>%
ungroup() %>%
rename(x_raw = x, y_raw = y) %>%
left_join(zeros) %>%
mutate(y_raw = ifelse(all_zero, 0, y_raw)) %>%
arrange(x_raw)
y_scale <- .25
dens_diag <- dens_df %>%
mutate(x1 = x_raw - y_scale * y_raw, x2 = x_raw + y_scale * y_raw,
y1 = x_raw + y_scale * y_raw, y2 = x_raw - y_scale * y_raw)
p <- ggplot(dens_diag) +
theme_ohara() +
geom_segment(aes(x = x1, y = y1, xend = x2, yend = y2),
color = 'grey80', alpha = 1, size = .1) +
geom_path(aes(x = x1, y = y1), color = 'grey20', alpha = 1, size = .1) +
geom_path(aes(x = x2, y = y2), color = 'grey20', alpha = 1, size = .1) +
geom_vline(xintercept = c(0, 1), color = 'grey20', size = .1) +
geom_hline(yintercept = c(0, 1), color = 'grey20', size = .1) +
facet_grid(tx_clean ~ str_clean) +
coord_fixed(xlim = c(0, 1), ylim = c(0, 1), expand = c(.02, 0)) +
geom_vline(data = mean_df, aes(xintercept = mean), color = 'red', size = .2) +
geom_hline(data = mean_df, aes(yintercept = mean), color = 'red', size = .2) +
geom_point(data = mean_df, aes(x = mean, y = mean), color = 'red') +
### add labels but inset just a bit
scale_x_continuous(breaks = c(0.1, 0.9), labels = c(0, 1)) +
scale_y_continuous(breaks = c(0.1, 0.9), labels = c(0, 1)) +
theme(axis.title = element_blank(),
legend.position = 'none',
panel.grid.major = element_blank(),
axis.ticks.length = unit(0, 'cm'),
strip.background = element_blank(),
strip.text.x = element_text(angle = 90, hjust = 0),
strip.text.y = element_text(angle = 0, hjust = 0))
ggsave('vuln_summary_plot_violin.png', width = 6, height = 6, dpi = 300)
knitr::include_graphics('vuln_summary_plot_violin.png')

taxa <- vuln_df$taxon %>% unique() %>% sort()
strs <- vuln_df$stressor %>% unique() %>% sort()
mean_df <- vuln_df %>%
group_by(tx_clean, str_clean) %>%
summarize(v_mean = mean(vuln, na.rm = TRUE),
v_sd = sqrt(mean(vuln_sd^2)),
n_spp = n_distinct(species),
.groups = 'drop') %>%
arrange(desc(v_mean)) %>%
mutate(tx_clean = fct_inorder(tx_clean) %>% fct_rev(),
str_clean = fct_inorder(str_clean))
x <- ggplot(mean_df, aes(x = str_clean, y = tx_clean)) +
theme_ohara() +
geom_tile(aes(fill = v_sd)) +
geom_point(color = 'white', shape = 15, size = 5) +
geom_point(aes(color = v_mean), shape = 15, size = 4) +
# geom_
scale_color_viridis_c(limits = c(0, 1)) +
scale_fill_gradient(low = 'white', high = 'grey20',
breaks = seq(0, .15, .03)) +
# scale_x_discrete(expand = c(0, 1.1)) +
coord_fixed() +
theme(axis.title = element_blank(),
# legend.position = 'none',
panel.grid.major = element_blank(),
axis.text.x = element_text(angle = 90, hjust = 1, vjust = .5)) +
labs(fill = 'std dev', color = 'mean',
title = 'Vulnerability by taxon and stressor')
ggsave('vuln_heatmap.png', width = 6, height = 4, dpi = 300)
knitr::include_graphics('vuln_heatmap.png')

str_summary_df <- vuln_df %>%
mutate(v_mean_all = mean(vuln, na.rm = TRUE)) %>%
group_by(str_clean) %>%
summarize(v_mean_all = first(v_mean_all),
v_mean = mean(vuln, na.rm = TRUE),
v_sd = sqrt(mean(vuln_sd^2)),
n_spp = n_distinct(species),
.groups = 'drop') %>%
arrange(v_mean) %>%
mutate(str_clean = fct_inorder(str_clean))
x <- ggplot(vuln_df %>%
mutate(str_clean = factor(str_clean, levels = str_summary_df$str_clean)),
aes(x = str_clean, y = vuln)) +
theme_ohara() +
geom_violin(scale = 'width', color = 'grey40', size = .2, fill = 'grey80', adjust = 2) +
geom_point(data = str_summary_df, aes(y = v_mean),
shape = 21, color = 'black', fill = 'red', size = 2) +
geom_hline(data = str_summary_df,
aes(yintercept = v_mean_all), color = 'red', size = .25) +
scale_y_continuous(expand = c(.01, .01)) +
coord_flip() +
theme(axis.title = element_blank(),
axis.text.x = element_text(angle = 90, hjust = 1, vjust = .5)) +
labs(title = 'Vulnerability distribution by stressor')
ggsave('vuln_dist_str.png', width = 4, height = 3, dpi = 300)
knitr::include_graphics('vuln_dist_str.png')

tx_summary_df <- vuln_df %>%
mutate(v_mean_all = mean(vuln, na.rm = TRUE)) %>%
group_by(tx_clean) %>%
summarize(v_mean_all = first(v_mean_all),
v_mean = mean(vuln, na.rm = TRUE),
v_sd = sqrt(mean(vuln_sd^2)),
n_spp = n_distinct(species),
.groups = 'drop') %>%
arrange(v_mean) %>%
mutate(tx_clean = fct_inorder(tx_clean))
x <- ggplot(vuln_df %>%
mutate(tx_clean = factor(tx_clean, levels = tx_summary_df$tx_clean)),
aes(x = tx_clean, y = vuln)) +
theme_ohara() +
geom_violin(scale = 'width', color = 'grey40', size = .2, fill = 'grey80', adjust = 2) +
geom_point(data = tx_summary_df, aes(y = v_mean),
shape = 21, color = 'black', fill = 'red', size = 2) +
geom_hline(data = tx_summary_df,
aes(yintercept = v_mean_all), color = 'red', size = .25) +
scale_y_continuous(expand = c(.01, .01)) +
coord_flip() +
theme(axis.title = element_blank(),
axis.text.x = element_text(angle = 90, hjust = 1, vjust = .5)) +
labs(title = 'Vulnerability distribution by taxon')
ggsave('vuln_dist_tx.png', width = 4, height = 2, dpi = 300)
knitr::include_graphics('vuln_dist_tx.png')

str_summary_df <- vuln_df %>%
mutate(v_mean_all = mean(vuln, na.rm = TRUE)) %>%
group_by(str_clean) %>%
summarize(v_mean_all = first(v_mean_all),
v_mean = mean(vuln, na.rm = TRUE),
v_sd = sqrt(mean(vuln_sd^2)),
n_spp = n_distinct(species),
.groups = 'drop') %>%
arrange(v_mean) %>%
mutate(str_clean = fct_inorder(str_clean))
x <- ggplot(vuln_df %>%
mutate(str_clean = factor(str_clean, levels = str_summary_df$str_clean)),
aes(x = str_clean, y = vuln)) +
theme_ohara() +
geom_boxplot(color = 'grey40', size = .2, fill = 'grey80', outlier.size = .2) +
geom_point(data = str_summary_df, aes(y = v_mean),
shape = 21, color = 'black', fill = 'red', size = 1.3) +
geom_hline(data = str_summary_df,
aes(yintercept = v_mean_all), color = 'red', size = .25) +
scale_y_continuous(expand = c(.01, .01)) +
coord_flip() +
theme(axis.title = element_blank(),
axis.text.x = element_text(angle = 90, hjust = 1, vjust = .5)) +
labs(title = 'Vulnerability distribution by stressor')
ggsave('vuln_dist_str_box.png', width = 4, height = 3, dpi = 300)
knitr::include_graphics('vuln_dist_str_box.png')

tx_summary_df <- vuln_df %>%
mutate(v_mean_all = mean(vuln, na.rm = TRUE)) %>%
group_by(tx_clean) %>%
summarize(v_mean_all = first(v_mean_all),
v_mean = mean(vuln, na.rm = TRUE),
v_sd = sqrt(mean(vuln_sd^2)),
n_spp = n_distinct(species),
.groups = 'drop') %>%
arrange(v_mean) %>%
mutate(tx_clean = fct_inorder(tx_clean))
x <- ggplot(vuln_df %>%
mutate(tx_clean = factor(tx_clean, levels = tx_summary_df$tx_clean)),
aes(x = tx_clean, y = vuln)) +
theme_ohara() +
geom_boxplot(color = 'grey40', size = .2, fill = 'grey80', outlier.size = .2) +
geom_point(data = tx_summary_df, aes(y = v_mean),
shape = 21, color = 'black', fill = 'red', size = 1.3) +
geom_hline(data = tx_summary_df,
aes(yintercept = v_mean_all), color = 'red', size = .25) +
scale_y_continuous(expand = c(.01, .01)) +
coord_flip() +
theme(axis.title = element_blank(),
axis.text.x = element_text(angle = 90, hjust = 1, vjust = .5)) +
labs(title = 'Vulnerability distribution by taxon')
ggsave('vuln_dist_tx_box.png', width = 4, height = 2, dpi = 300)
knitr::include_graphics('vuln_dist_tx_box.png')

tx_summary_df <- vuln_df %>%
group_by(tx_clean) %>%
mutate(v_mean_all = mean(vuln, na.rm = TRUE)) %>%
group_by(tx_clean, str_clean) %>%
summarize(v_mean_all = first(v_mean_all),
v_mean = mean(vuln, na.rm = TRUE),
v_sd = sqrt(mean(vuln_sd^2)),
n_spp = n_distinct(species),
.groups = 'drop') %>%
arrange(desc(v_mean)) %>%
mutate(tx_clean = fct_inorder(tx_clean),
str_clean = fct_inorder(str_clean) %>% fct_rev())
plot_df <- vuln_df %>%
mutate(tx_clean = factor(tx_clean, levels = levels(tx_summary_df$tx_clean)),
str_clean = factor(str_clean, levels = levels(tx_summary_df$str_clean))) %>%
# slice_sample(prop = .5) %>%
rowwise() %>%
mutate(v = list(rnorm(n = 30, mean = vuln, sd = vuln_sd))) %>%
unnest(v) %>%
mutate(v = case_when(v < 0 ~ 0, v > 1 ~ 1, TRUE ~ v)) %>%
ungroup()
x <- ggplot(plot_df, aes(x = str_clean, y = v)) +
theme_ohara(base_size = 7) +
geom_boxplot(color = 'grey40', size = .2, fill = 'grey80',
outlier.size = .05, outlier.color = 'grey70') +
# geom_violin(scale = 'width', color = 'grey40', size = .2, fill = 'grey80') +
# geom_jitter(color = 'black', size = .2, alpha = .1, width = .4, height = 0) +
geom_hline(data = tx_summary_df,
aes(yintercept = v_mean_all), color = 'red', size = .25) +
geom_point(data = tx_summary_df, aes(y = v_mean),
shape = 21, color = 'black', fill = 'red', size = 1) +
scale_y_continuous(expand = c(.01, .01)) +
coord_flip() +
theme(# strip.background = element_blank(),
strip.text = element_text(face = 'bold', hjust = 0),
axis.title.y = element_blank(),
axis.text.x = element_text(angle = 90, hjust = 1, vjust = .5)) +
labs(# title = 'Vulnerability distribution by taxon and stressor',
y = 'Vulnerability') +
facet_wrap( ~ tx_clean, nrow = 3)
ggsave('vuln_dist_tx_3x4panel.png', width = 6, height = 6, dpi = 300)
knitr::include_graphics('vuln_dist_tx_3x4panel.png')

# tx_summary_df <- vuln_df %>%
# group_by(tx_clean) %>%
# mutate(v_mean_all = mean(vuln, na.rm = TRUE)) %>%
# group_by(tx_clean, str_clean) %>%
# summarize(v_mean_all = first(v_mean_all),
# v_mean = mean(vuln, na.rm = TRUE),
# v_sd = sqrt(mean(vuln_sd^2)),
# n_spp = n_distinct(species),
# .groups = 'drop') %>%
# arrange(desc(v_mean)) %>%
# mutate(tx_clean = fct_inorder(tx_clean),
# str_clean = fct_inorder(str_clean) %>% fct_rev())
#
# plot_df <- vuln_df %>%
# mutate(tx_clean = factor(tx_clean, levels = levels(tx_summary_df$tx_clean)),
# str_clean = factor(str_clean, levels = levels(tx_summary_df$str_clean))) %>%
# # slice_sample(prop = .05) %>%
# rowwise() %>%
# mutate(v = rnorm(n = 1, mean = vuln, sd = vuln_sd),
# v = case_when(v < 0 ~ 0, v > 1 ~ 1, TRUE ~ v)) %>%
# ungroup()
x <- ggplot(plot_df %>% mutate(str_clean = fct_rev(str_clean)),
aes(x = str_clean, y = v)) +
theme_ohara(base_size = 7) +
### manually create panel background and gridlines
scale_x_discrete() +
annotate('rect', xmin = .5, xmax = 22.5, ymin = 0, ymax = 1,
color = NA, fill = 'grey95') +
geom_hline(yintercept = c(.25, .5, .75, 1), color = 'white', size = .2) +
geom_vline(xintercept = 1:22, color = 'white', size = .2) +
geom_hline(yintercept = 0, color = 'grey20', size = .2) +
scale_y_continuous(breaks = c(0.05, 0.50, 0.95), labels = c('0.0', '0.5', '1.0')) +
### boxplot
geom_boxplot(color = 'grey40', size = .2, fill = 'grey80', outlier.size = .05, outlier.color = 'grey70') +
geom_hline(data = tx_summary_df,
aes(yintercept = v_mean_all), color = 'red', size = .1) +
geom_point(data = tx_summary_df, aes(y = v_mean),
shape = 21, color = 'black', fill = 'red', size = 1) +
theme(strip.background = element_blank(),
strip.text.y = element_text(face = 'bold', hjust = 0, angle = 0),
axis.title.x = element_blank(),
axis.text.x = element_text(angle = 90, hjust = 1, vjust = .5)) +
labs(# title = 'Vulnerability distribution by taxon and stressor',
y = 'Vulnerability') +
facet_wrap( ~ tx_clean, ncol = 1, strip.position = 'right')
ggsave('vuln_dist_tx_12x1panel.png', width = 4, height = 8, dpi = 300)
knitr::include_graphics('vuln_dist_tx_12x1panel.png')
